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Abstract 

The ground-state structures and ferroelectric properties of Ndat(./V=2-52) have been investigated 
by a combination of density-functional theory (DFT) in the generalized gradient approximation 
(GGA) and an unbiased global search with the guided simulated annealing. It is found that the 
electric dipole moment (EDM) exists in the most of Nb^v and varies considerably with their sizes. 
And the larger Nbjv {N >25) prefer the amorphous packing. Most importantly, our numerical EDM 
values of Nbjy (N >38) exhibit an extraordinary even-odd oscillation, which is well consistent with 
the experimental observation, showing a close relationship with the geometrical structures of NbAr. 
Finally, an inverse coordination number (ICN) function is proposed to account for the structural 
relation of the EDM values, especially their even-odd oscillations starting from Nb38. 

PACS numbers: 36.40.Mr, 31.15.Ar, 36.90.+f, 73.22.-f 
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I. INTRODUCTION 



Transition metal clusters have long been of considerable interest due to their importance 
in fundamental researches and tremendous potential technical applications, among which 
niobium cluster is one of the most thoroughly studied in both of their chemical and physical 
properties. Even so, the fully understanding of their dramatic size-dependent properties is 
still a great challenge. Recently, the ferroelectricity (FE) in free niobium clusters (Nbjv, 
iV= 2-150) has been experimentally found at low temperatures 1,2 , showing the existence of 
electric dipole moments (EDMs), which, more importantly, exhibit a pronounced even-odd 
oscillation for Nbjv (N >38). Such an interesting phenomenon has attracted immediate at- 
tention because the FE has never been found in single-element bulk materials and neither in 
metals. For example, a purely electronic mechanism has been proposed, considering the FE 
of Nb^Y to be caused by the electron correlations 3 . However, the large size-dependent EDMs, 
especially the remarkable even-odd oscillation starting from N = 38, strongly support ex- 
istence of an intimate relation between the FE and structures of Nb^. Since the existing 
techniques cannot conclusively determine the geometrical structures of clusters, the numer- 
ical simulation on them becomes a useful method to provide more valuable information on 
their structures and novel physical properties. 

Several first-principle calculations on Nb^ up to 23 atoms had been performed 4-6 . It was 
concluded that the icosahedral growth is not favored for Nbjv 6 . Other studies made on the 
smaller Nbjv (N <10) found the similar lowest-energy structures with the high coordinated 
configurations. Recently, a close relationship between the asymmetrical geometrical struc- 
tures of Nbjv (N <15) and their EDMs has been revealed through a first-principle study 7 . 
However, there has been no DFT study on the larger Nb^ and their ferroelectric properties. 

Therefore, in this paper, we will address the growth pattern of Nbjy {N <52) and its 
close relation with cluster FE. The relevant structures of Nb^ are obtained by an empirical 
global optimization combined with the DFT relaxation, from which the EDMs are then 
calculated by the DFT. The details of the calculated methods are described in Sec. II. The 
structural and ferroelectric properties of Nbjv are presented in Sec. III. Section IV contains 
the concluding remarks of our work. 
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II. COMPUTATIONAL METHOD 



Instead of the presumed symmetry constraints, an unbiased global search on the potential 
energy surface of clusters is preformed by the guided simulated annealing (GSA) method 8 ' 9 , 
which incorporates a guiding function (GF) in the traditional simulated annealing to find the 
global energy minimum. In our calculations, the density of atoms and the second moment 
of the mass distribution 10 are adopted as the GFs. The interaction between Nb atoms is 
represented by an empirical many-body potential 11 ' 12 , which has been successfully applied 
to predict the equilibrium geometries of Nbjv {N <14) 13 . A number of lower-energy isomers 
are thereby generated, representing different local stable states in the phase space. 

Then, the DMol3 package 14 based on DFT is used to further optimize the cluster struc- 
tures by selecting at least five isomers with different symmetries for a given size from the 
above step, in which a relativistic effective core potential (RECP) 15 ' 16 and a double nu- 
merical basis including a d-polarization function are chosen to do the electronic structure 
calculations. The RECP was generated by fitting all-electron Hartree-Fock results, which 
has been included in the DMol3 package as one of the powerful method treating the core 
electrons. The density functional is treated in GGA with spin polarization, and Perdew and 
Wang exchange-correlation function is used 17 . Geometrical structure optimization is per- 
formed with the Broyden-Fletcher-Goldfarb-Shanno algorithm 18 . A convergence criterion of 
10~ 3 a.u. on the gradient and displacement, and 1CT 5 a.u. on the total energy is used in 
the optimization. The accuracy of the current computational scheme has been checked by 
benchmark calculations on the Nb atom and its bulk solid 19 . In addition, a comparison has 
been made on the results of Nb^ (A=2-19, 43-45) obtained by both RECP and all-electron 
calculations with scalar relativistic corrections (AESRC), showing a reasonable consistency 
between them in predicting the structural order and corresponding properties of Nb^ clus- 
ters. With this strategy, we have optimized the equilibrium structures of Nb^r up to 52 
atoms in a reliable and efficient way, though we cannot strictly rule out other energetically 
more favorable structures. 

III. RESULTS 

The calculated binding energies (BEs) per atom are shown in Fig. 1, in which the isomers 
and the ground-state structures of some clusters are also included. For small Nb N (N=2- 
12) 20 , our lowest-energy configurations and isomeric order are similar to those found in 
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Ref. 6 except for Nbs, which has several competitive isomers, such as a capped distorted 
pentagonal bipyramid and a C 2t , capped octahedron, within an energy interval of 0.1 eV. 
Though the atomic packing of the small clusters is gradually changed from one size to 
another, there exist several structural transitions for medium-sized clusters of N— 13-24, 
e.g. from the icosahedron to the icositetrahedron, as illustrated in Fig. 1. Nt>i3 is a 
distorted icosahedron. By adding one atom to Nbi3, one planar pentagon in it will be 
changed into the first hexagonal ring, forming the structure of Nbu. Again, one more atom 
makes the opposite planar pentagon in Nb i4 to become the second hexagonal ring, leading 
to the equilibrium icositetrahedral structure of Nbi 5 , which is about 0.26 eV and 0.89 eV 
lower in energy than the distorted bulk bcc type structures with approximate Oh and D 2 h 
symmetries, respectively (see Fig. 1). The structures of Nbi6-Nbig are obtained by adding 
1, 2, and 3 bottom atoms to Nbis, holding still the two hexagonal rings (distorted). Nbig 
is a double icosahedral structure with two core atoms. From Nb 2 o to Nb 22 , the structures 
evolve into a slightly distorted double-interpenetrating icositetrahedron with three parallel 
hexagonal rings, which are obtained by adding one, two, and three atoms to the three 
pentagons of Nbig, respectively. Nb 2 3 is formed by adding one atom to the bottom of Nb 22 , 
and by further adding an edge-capped atom, Nb 2 4 is obtained. 

Now, we pay more attention to the larger Nbjv (A=25-52) by the same computational 
scheme. The obtained results suggest that they prefer the amorphous packing, making the 
compact oblate-spherical configurations dominant. The average coordination number and 
bond length as well as the ratio of body atom number to surface atom one (N b /Nf) are 
shown in Fig. 2 for 25< N <52. Here, the nearest-neighbor bond length is truncated at 3.3 
A, which is obtained by calculations of atom pair-correlation function g(r) for all clusters. 
It can be seen clearly from Fig. 2 that both of the average coordination number and bond 
length increase non-monotonically with size, showing a close similarity to each other. The 
apparent size- dependent geometrical structures of Nbjv may play an important role in their 
physical properties. For example, the size- dependent deviation of the bond length from the 
bulk value may induce the symmetry breaking and so the emergence of EDM. It is also 
worth noticing that the dramatic even-odd oscillation behavior appears in Nb/Nf of Nb^ 
(N >40), suggesting the growth pattern may be dualistic, which is different from even to 
odd Nbjv, leading to a probable even-odd variation of their physical properties. 

Based upon the above structures, the second order difference A 2 E(A) of BE varying 
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with size is calculated and shown in the inset of Fig. 1, which evidently displays a good 
agreement between both results from RECP and AESRC. The variation of A 2 E(iV) shows 
that the magic number should emerge at N=7, 13, 15, 17, 22, 24, 27, 29, 31, 33, 35, 37, 41, 
44, 47, and 50, which agrees well with the peaks in the abundant spectra of Nb^r at N—7, 
13, 15, 22, 29, and 33 21 . 

The extraordinary FE of Nbjv observed in the experiment 1 is particular interesting. The 
calculated EDM curves are shown in Fig. 3 (a), from which we find again that the change of 
the RECP results agrees with that of the AESRC ones. In the most cases, the RECP EDM 
values are smaller than the AESRC ones. However, the positions of local maximum and 
minimum EDMs obtained by RECP and AESRC are in a good agreement with each other, 
which are also consistent with the experimental observations. For example, we obtained 
theoretically the local maximum EDMs at iV=14, 18, 20, 24, and local minimum ones at 
iV=13, 15, 19, 22 for the medium-size Nbjv (iV = 13-24). More importantly, a markedly 
even-odd EDM oscillation has been also reproduced theoretically, starting from Nb 38 . 

Obviously, the EDM of a niobium cluster is determined by its asymmetrical charge dis- 
tribution (CD), which should be greatly influenced by its geometrical structure, such as 
its shape, surface atom number, and inter-atomic distances deviated from the bulk values, 
all of which induce a deformed CD. For example, the lowest-energy geometrical structure 
of Nb 6 is a distorted prism, which is more stable than the octahedron only by 0.03 eV. 
The approximate Oh symmetry of the latter prohibits appearance of the EDM, while the 
former has an EDM of 0.2721 D due to inverse symmetry breaking. That is to say, the 
different isomers have different EDMs even if they are energetically close, showing also the 
geometrical structure of a niobium cluster has a very important effect on its EDM. With 
increase of atoms, there exist some metastable isomers close to the lowest-energy structure 
in the most of Nb^, whose EDM amplitudes may have a larger fluctuation. In order to 
detect further the isomer's effect on the EDM, we have given the results of some larger Nb^ 
in Table 1, showing clearly the EDMs of different isomers fluctuate greatly. Typically, the 
variation of EDM between isomers, within a vibrational temperature of 100 K 22 , at a given 
size is about 0.3 - 0.9 D, which is obviously larger than the even-odd oscillation amplitude of 
EDMs between neighboring Nbjv (about 0.1-0.3 D). For example, the ground-state structure 
of Nb 5 i and its first close-lying isomer are close to each other in energy, but their EDMs 
differ heavily from each other, varying from 0.3020 D to 1.1736 D. As shown in Ref. 1, the 
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electric deflection of Nb^ has been measured at the finite temperatures, showing the clearer 
even-odd oscillation starting from Nb 38 at a lower temperature of 20 K. With temperature 
increasing, the oscillation amplitude becomes smaller and smaller, and finally disappears 
above 100 K. In our calculations, the isomer's influence on the EDM value plays also a "fi- 
nite temperature" effect. The better the structure optimization is, the clearer the even-odd 
oscillation behavior. In addition, we also found from Table 1 that the average coordination 
number of the ground-state structures is the largest in the isomers we considered, showing 
that the ground-state structures of Nb^ prefer the compact configurations. 

In order to further find what reason to cause the anomalous EDMs of Nb^, we have 
introduced a concept of the effective charge to characterize the bonding variations in a 
cluster, which could be obtained by calculation of the coordination number because it can 
reflect, to some extent (although not very precisely), the size- dependent structures of Nbjv. 
Due to the un-equivalent geometrical surroundings of each atom in a cluster, the electron 
density should not distribute uniformly between two nearby atoms forming a chemical bond, 
leading to different effective charges on the atoms of a Nb^ and so emergence of its EDM. 
It is obvious that the effective charge on an atom decreases with increase of its coordination 
number, i.e. an atom with higher coordination number should possess less effective charge. 
Therefore, we adopt simply an inverse coordination number (ICN) as a weight index to 
quantify the effective CD in a Nb^, which could be represented by a function F(N), defined 
as: 



F(N) = 



N 



B £ 



Ri / Zi + Rj / Zj 



1/Zi + 1/Zj 



hJ = 1 

r ij<R cu t 



1 N 

where B = ^ Zi is the total bonding number with Zi the coordination number of the ith 

i=i 

atom. R C ut is the cutoff distance (3.3 A), and Ri is the position vector of the ith atom (the 
coordinate origin is set at the mass center of the cluster). 

The F(N) values of Nb^ (N— 2-52) are shown in Fig. 3 (b), which displays almost 
the same variation behavior as the EDMs, demonstrating it correctly describes the CD 
deformation in Nb^. For example, for the smaller Nb^, the ICN function reproduces again 
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the local maxima and minima of EDMs at N=6, 11, 18, 20, 24, 28, 30 and at iV = 4, 7, 10, 
13, 15, 17, 19, 29, respectively. Specially, for the larger Nb^ with N >38, the ICN values 
are enhanced for even clusters, but suppressed for odd ones, reproducing the extraordinary 
even-odd EDM oscillation shown in Fig. 3(a), which clearly indicates that the ICN function 
indeed can be used to characterize the effective CD in Nb^ qualitatively although it is so 
simple and elegant. 

The close relation between the Nb^ structures and their corresponding EDMs can also be 
identified by visualization of the spatially deformed CD in the clusters, which is defined as the 
total charge density minus the density of the isolated atoms. Thus, the regions with positive 
deformation charge density will indicate formation of the bonds, while the negative regions 
indicate electron loss. For example, the bonding characters of Nbig and Nb2o with isodensity 
surface of value 0.05 e/au 3 are illustrated in Fig. 4, which clearly shows the difference 
between their deformation charge densities. A slightly distorted double-icosahedral structure 
of Nbig induces its charge isocontour with an approximate D 5 h symmetry, showing an almost 
zero EDM. However, the obviously different CD from top to bottom of Nb 20 shown in Fig. 
4(b) gives rise to its rather larger EDM, denoted by a yellow arrow. So, we conclude at this 
point that the different structures of Nb^ will induce much different spatial CD, leading to 
different EDM values of Nb^v- 

IV. CONCLUSION 

In summary, the equilibrium geometries, relative stabilities, and EDMs of Nb^ (N = 
2-52) have been calculated by a combination of empirical interaction model and the DFT 
optimization. The more special attention has been paid on the effects of the Nbjv structures 
on their EDMs. It is found that no one Nb^ mimics the bulk structure and the compact 
oblate-spherical amorphous structures are preferable for the larger Nb^r (N >25). Inter- 
estingly, the size-dependent structures of Nb^r are found to be an intrinsic origin to induce 
their unordinary FE. Our study shows that the EDM does exist in the most of niobium 
clusters and has a close relationship with their geometrical structures. A simple ICN func- 
tion is proposed to account for the anomalous size- and structure-dependent variations in 
the EDMs of Nb^. A good agreement between the ICN function and the theoretical values 
of EDMs as well as the experiment demonstrates the geometrical structure of Nb^r has an 
important effect on its ferroelectric property. 
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For a niobium atom, the ionization potential and electron affinity obtained by RECP are 7.011 
eV and 0.758 eV, respectively, which agree well with the experimental results: 6.76 eV and 0.89 
eV. The calculated lattice parameter of 3.30 A and cohesive energy of 7.26 eV per atom are 
also consistent with the experimental data of 3.30 A and 7.57 eV/atom, respectively. 
The BE and bond length of Nb 2 , obtained by RECP (AESRC), are 3.38 (4.88) eV and 2.18 
(2.15) A, respectively. The BE of AESRC agrees well with the experimental value of 4.86±0.02 
eV [M. D. Morse, Chem. Rev. 86, 1049 (1986)], but that of RECP underestimates it. However, 
as we shall show, the structural and energetic features for Nb^ obtained by the two methods 
are in a satisfactory agreement. 

M. Sakurai, K. Watanabe, K. Sumiyama, and K. Suzuki, J. Chem. Phys. Ill, 235 (1999). 
The vibrational temperature is calculated in terms of the energy difference by T=2AE/(3N- 
G)ks, where ks is Boltzmann constant. 
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TABLE I: Relative energy, the EDM, and the average coordination number (Z) of the ground-state 
and its first two isomers of Ndat at iV=38-40, 50-52. The energetic differences between the isomers 
are at the limits of the accuracy of our DFT-GGA approach. 



N 


E (eV) 


EDM (D) 


Z 


N 


E (eV) 


EDM (D) 


Z 


38a 





0.5290 


8.158 


50a 





0.4677 


8.360 


38b 


0.074 


0.3511 


7.894 


50b 


0.203 


0.4427 


7.880 


38c 


0.268 


0.2675 


7.947 


50c 


0.452 


0.5145 


7.840 


39a 





0.2018 


8.231 


51a 





0.3040 


8.275 


39b 


0.072 


0.4917 


7.846 


51b 


0.025 


1.1736 


8.196 


39c 


0.184 


0.5935 


7.897 


51c 


0.163 


1.1473 


8.196 


40a 





0.5758 


8.050 


52a 





0.5696 


8.269 


40b 


0.176 


0.2920 


8.025 


52b 


0.142 


0.8694 


8.038 


40c 


0.257 


0.5215 


7.700 


52c 


0.325 


0.3374 


7.750 
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Figure Captions 



FIG. 1. (Color online) The binding energy (BE) per atom of Nb^ vs cluster size. The inset is 
the size-dependent second order difference of BE. The RECP and AESRC results are shown 
by circles and triangles, respectively. The vertical dashed lines represent the positions at N 
= 6, 13, 15, 19, and 22. Geometrical structures of the different isomers for Nb 6 and Nbi 5 are 
also shown in the order of their BE values from top to bottom. The structural transitions 
are illustrated by the ground-state structures of Nbi3, Nbis, Nbig, and Nb 2 2- 

FIG. 2. (Color online) Average coordination number (average bond length) vs cluster size, 
guided for eyes with solid and dashed lines, respectively. The inset shows the ratio of body 
atom number to surface atom one, denoted as N b /Nf. 

FIG. 3. (Color online) (a) DFT values of EDMs vs cluster size, calculated by two different 
methods, i.e. RECP (circles) and AESRC (triangle), (b) Dependence of the ICN function 
on cluster size. For clarity, the values of even and odd clusters are shown by open and filled 
circles, respectively. 

FIG. 4. (Color online) The deformation densities for (a) Nbig, (b) Nb 20 , in which the direction 
of the EDM is denoted by the yellow arrow. The isodensity surface corresponds to 0.05 
e/a.u. 3 . 



11 



